# params <- c('sigy','Na','N1',
#             'alpha1','alphaa','alphar','alphal',
#             'beta1', 'betaa','betar', 'betal',
#             'r')


# G and P only for some initial runs to check whether they make sese...
params <- c('sigy','Na',
            'alpha1','alphaa','alphar','alphal',
            'beta1', 'betaa','betar', 'betal')
            # 'G','P')

# inits <- function(){list(tauy = 1, 
#                           Na = c(1000, 1000, 1092.23, 1100.01, 1234.32, 1460.85, 1570.38, 1819.79,
#                           1391.27, 1507.60, 1541.44, 1631.21, 1628.60, 1609.33, 1801.68, 1809.08, 1754.74,
#                           1779.48, 1699.13, 1681.39, 1610.46, 1918.45, 1717.07, 1415.69, 1229.02, 1082.02,
#                           1096.61, 1045.84, 1137.03, 981.1, 647.67, 992.65, 968.62, 926.83, 952.96, 865.64),
#                           N1 = rep(400,36),
#                           alpha1 = 1, alphaa = 2, alphar = -2, alphal = -4, 
#                           beta1 =-2, betaa = 0.1, betar = -0.7, betal = -0.3)}

# round() for discrete distribution
# inits <- function(){list(tauy = 1,
#                          Na = round(c(1000, 1000, 1092.23, 1100.01, 1234.32, 1460.85, 1570.38, 1819.79,
#                                 1391.27, 1507.60, 1541.44, 1631.21, 1628.60, 1609.33, 1801.68, 1809.08, 1754.74,
#                                 1779.48, 1699.13, 1681.39, 1610.46, 1918.45, 1717.07, 1415.69, 1229.02, 1082.02,
#                                 1096.61, 1045.84, 1137.03, 981.1, 647.67, 992.65, 968.62, 926.83, 952.96, 865.64))/sc,
#                          alpha1 = 1, alphaa = 2, alphar = -2, alphal = -4,
#                          beta1 =-2, betaa = 0.1, betar = -0.7, betal = -0.3)}

# inits <- function(){list(tauy = 1, 
#                          Na = round(c(1000, 1000, 1092.23, 1100.01, 1234.32, 1460.85, 1570.38, 1819.79,
#                                       1391.27, 1507.60, 1541.44, 1631.21, 1628.60, 1609.33, 1801.68, 1809.08, 1754.74,
#                                       1779.48, 1699.13, 1681.39, 1610.46, 1918.45, 1717.07, 1415.69, 1229.02, 1082.02,
#                                       1096.61, 1045.84, 1137.03, 981.1, 647.67, 992.65, 968.62, 926.83, 952.96, 865.64))/sc,
#                          alpha1 = 0.54, alphaa = 2, alphar = -2, alphal = -4, 
#                          beta1 = -0.19, betaa = 0.1, betar = -0.7, betal = -0.3)}
 
inits <- function(){list(tauy = 1,
                         Na = round(c(1000, 1000, 1092.23, 1100.01, 1234.32, 1460.85, 1570.38, 1819.79,
                                      1391.27, 1507.60, 1541.44, 1631.21, 1628.60, 1609.33, 1801.68, 1809.08, 1754.74,
                                      1779.48, 1699.13, 1681.39, 1610.46, 1918.45, 1717.07, 1415.69, 1229.02, 1082.02,
                                      1096.61, 1045.84, 1137.03, 981.1, 647.67, 992.65, 968.62, 926.83, 952.96, 865.64)/sc),
                         alpha1 = 1, alphaa = 2, alphar = -1.106, alphal = -4,
                         beta1 = -0.19, betaa = 0.1, betar = -0.299 , betal = -0.3)}

 
### FULL DA posterior means
# alpha1     0.5540 6.888e-02 
# alphaa     1.5678 6.338e-02 
# alphal    -4.5760 3.534e-02 
# alphar    -1.1723 6.579e-02 
# beta1     -0.1913 5.668e-02 
# betaa     -0.2465 3.786e-02 
# betal     -0.3650 3.981e-02 
# betar     -0.3360 3.254e-02 
# sigy   30755.8259 8.700e+03 

# Na[3]    993.9475 26.52
# Na[13]  1789.2550 52.95
# Na[23]  1469.4150 49.46
# Na[33]   975.9720 59.51
Na = round(c(1000, 1000, 1092.23, 1100.01, 1234.32, 1460.85, 1570.38, 1819.79,
             1391.27, 1507.60, 1541.44, 1631.21, 1628.60, 1609.33, 1801.68, 1809.08, 1754.74,
             1779.48, 1699.13, 1681.39, 1610.46, 1918.45, 1717.07, 1415.69, 1229.02, 1082.02,
             1096.61, 1045.84, 1137.03, 981.1, 647.67, 992.65, 968.62, 926.83, 952.96, 865.64)/sc)